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Abstract.  Given  a planar  triangulation,  the  goal  is  to  select  elevations  at  its  vertices  so  that  result- 
ing piecewise-linear  triangulated  surface  approximates  specified  elevations  using  the  Li-norm  as  the 
primary  measure-of-fit.  Several  suboptimal  algorithms  relating  to  that  problem  have  been  devised 
and  implemented,  and  are  described  here,  as  part  of  TIN  prototype  software  for  terrain  modeling 
in  the  context  of  distributed  simulation  as  well  as  elevation  data  editing.  L\ -approximation  acts 
as  a median  filter  and  thus  provides  advantages  for  “bare  earth”  representation.  For  added  flexi- 
bility, a “sign-oriented  scaling”  scheme  has  been  developed,  which  generalizes  standard  norm-based 
measures-of-fit.  L\ -approximation  is  either  tackled  directly  or  by  iterating  suitably  weighted  L 2- 
approximations.  Numerical  experiments  are  reported. 

Key  words,  bare  earth,  distributed  simulation,  Li-approximat,ion,  ./^-approximation,  linear  pro- 
gramming, median  filter,  terrain  surface  modeling,  TIN,  triangulation 


1.  Introduction 

“Triangulated  irregular  networks  (TINs)”  are  playing  an  increasingly  prominent  role  in  mod- 
eling terrain  surfaces  for  purposes  such  as  animation,  distributed  simulation,  and  geodetic 
volumetries  (e.g.  [3]).  Typically,  there  is  a rectangular  map  area  together  with  a set  R of 
location  points  p — (xp,yp)  on  that  map  each  with  elevation  zp.  Terrain  is  to  be  modeled  as 
a surface  based  on  those  elevation  data. 

A triangulated,  irregular  network  is  a planar  triangulation  which  arises  as  the  footprint  of  a 
triangulated  surface  z = z(xpy ),  that  is,  a piecewise  linear  surface  consisting  of  flat  triangles 
joined  continuously  together  in  space.  That  surface  is  referred  to  as  TIN -surface.  Under 
triangulation , we  understand  any  covering  of  the  map  by  triangles  such  that  the  triangle 
interiors  do  not  overlap  and  any  two  triangle  edges  which  meet  in  more  than  one  point  are 
identical. 

Typically,  the  vertices  of  the  triangulation  have  previously  been  selected  as  critical  points 
based  on  “greedy”  criteria.  We  denote  the  set  of  those  vertices  by  S.  The  TIN  surface 
z = z(x,y)  is  uniquely  determined  by  the  elevations  zs  at  the  critical  points  s = ( xs,ys ) E S 
and  the  footprint  triangulation,  the  TIN'.  Indeed,  the  elevations  at  the  three  corners  of  any 
triangle  in  the  TIN  uniquely  determine  a plane  through  the  corresponding  elevation  points 
and,  therefore,  a triangle  in  space. 

The  object  then  is  to  approximate  given  elevation  data  by  a TIN  surface  - as  accurately  as 
possible  with  the  number  of  triangles  in  the  TIN  not  exceeding  a prescribed  upper  limit,  the 
“budget”.  Such  restrictions  are  necessary,  for  instance,  to  meet  the  high  speed  requirements 
of  real-time  visualization. 
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Generally,  there  are  two  aspects  of  this  optimization.  The  first  is  the  selection  of  the  TIN 
itself,  that  is.  of  the  critical  points  and  their  triangulation.  The  second  is  to  assign  elevations  at 
those  critical  points  so  as  to  optimize  some  measure-of-fit  to  the  given  elevations  while  keeping 
the  triangulation  fixed.  This  work  is  concerned  with  that  second  aspect  of  optimization. 

There  are  different  measures-of-fit  in  general  use  to  assess  the  accuracy  of  such  an  ap- 
proximation. Most  are  based  on  the  notion  of  a residual  at  each  given  elevation  location 
V = (iP,  yP ) e R , 

Tp  — Zp  Zp  z{ Xp , IJp ) , 

where  zp  denotes  the  elevation  z(xp,yp)  which  the  TIN  surface  assumes  at  location  {xp,yp). 
Two  commonly  considered  measures-of-fit  are  the  maximum  deviation 

MAX  = maxpefl|rp| 


and  the  root  mean  square 


RMS 


peR 


where  \R\  is  the  number  of  elevation  locations.  Both  quantities  relate  to  norms  of  vectors. 
MAX  is  an  instance  of  the 

Loo  ~ norm , 


also  called  Chebychev  norm.  RMS  is  a scaled  version  of  the 


L2  — norm 


or  least-squares  norm. 

There  is,  however,  a problem  with  L2-  and,  in  particular,  with  Loo-norms,  in  that  outliers 
among  the  elevation  data  caused  by  terrain  artifacts  such  as  rocks,  single  trees  and  bushes, 
as  well  as  processing  artifacts,  distort  the  shape  of  the  TIN  surface  - make  it  conform  less  to 
bare  earth  during  the  process  of  selecting  elevations  at  critical  points.  In  some  approaches, 
the  selection  of  a norm  may  also  affect  the  construction  of  the  TIN  itself. 

Those  drawbacks  of  the  commonly  considered  norms  have  led  the  authors  to  consider  the 

Li  — norm 

as  a measure-of-fit  for  terrain  surface  approximation  because  it  is  more  robust  with  respect 
to  artifacts,  and  may  even  be  used  for  their  detection  (see  [4]).  The  Li-norm  is  defined  as 
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the  sum  of  the  absolute  values  of  the  residuals.  One  might  actually  prefer  the  mean  absolute 
deviations 


MAD 


1 

W\ 


E K 

p£R 


Note  that  the  elevation  2:  = z(x,y)  of  the  TIN  surface  at  any  map  location  {x,y)  is  a 
linear  function  of  the  critical  elevations  zs,  s G S.  Consequently,  the  residuals  rp,  p £ R 
are  also  linear  functions  of  the  critical  elevations.  For  this  reason,  the  Li-approximation 
problem  is  a linear  programming  problem,  and  the  /^-approximation  problem  is  a linear  least 
squares  problem.  The  connection  between  Li-approximation  and  linear  programming  is  well 
established  (e.g.  [1] ) . 

More  precisely,  2 = z[x , y)  is  a linear  combination  of  three  elevations 


zs\i  zs 21  zs3i 

which  belong  to  vertices  si,  S2,  S3  of  a TIN  triangle  containing  location  ( x,y ).  The  coefficients 
of  that  linear  combination, 

zix,y)  — X\ZSl  + ^2As2  + A32S3, 

are  the  barycentric  coordinates  of  location  (x,  y)  with  respect  to  s*  = (xSi,ySi),  i — 1,  2,  3 and 
are  defined  by 


x = X1xSl  + X2xS2  + X3xS3 

V — Aiysx  + X2yS2  + ^3  yS3 

1 = Ai  + X2  + A3. 

The  barycentric  coordinates  assume  values  between  0 and  1 since  the  location  (x,  y)  lies  in 
the  triangle  in  question.  The  barycentric  coordinates  play  a major  role  in  approximation  by 
triangulated  surfaces. 

While  linear  programming  problems  are  well  understood  and  excellent  solution  methods 
exist  even  for  large  numbers  of  variables,  the  Li-approximation  problem  considered  here  typ- 
ically exceeds  the  size  that  can  be  handled  by  those  methods.  Suboptimal  solution  methods 
were,  therefore,  developed.  Straight  /^-approximation,  however,  can  be  solved  fully  optimally 
on  an  iterative  basis  described  below. 

One  method  of  elevation  optimization  with  respect  to  a specific  measure-of-fit,  is 

vertex  — based  elevation  optimization. 

With  any  critical  point  - vertex  in  the  TIN  we  associate  the  union  of  all  triangles  in  the 
triangulation  which  contain  that  vertex,  calling  it  the  star  of  that  critical  point.  The  star  is  a 
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simply-connected  polygon  bounded  by  the  edges  between  adjacent  “neighbor”  vertices  of  the 
critical  point  that  centers  the  star. 

We  will  require  the  triangulation  to  be 

star  — complete , 

that  is,  every  star  in  the  triangulation  contains  at  least  one  elevation  location  (. xp,yp ),  p G R 
in  its  interior.  Without  this  condition,  the  elevation  at  the  center  of  an  empty  star  could  be 
chosen  arbitrarily  without  affecting  any  of  the  norm-based  measures-of-fit  considered  here.  In 
that  sense,  the  approximation  problem  would  be  ill-defined. 

A single  step  of  a vertex-based  iteration  then  consists  of  adjusting  the  elevation  at  the 
center  of  a star  so  as  to  minimize  the  norm  of  choice  of  the  residuals  in  the  interior  of  the 
star  while  keeping  the  elevations  at  all  other  vertices  fixed.  Each  vertex  elevation  is  adjusted 
in  turn  in  a pass  through  all  critical  points.  Such  a pass  is  then  repeated  iteratively  until  the 
resulting  improvements  become  acceptably  small. 

For  Li-  and  L2-norms,  the  basic  vertex-based  optimization  step  can  be  executed  readily 
in  closed  form.  In  the  L2  case,  a convex  quadratic  of  one  variable  has  to  be  minimized.  In  the 
Li  case,  minimization  of  a piecewise  linear  convex  function  of  one  variable  is  required.  This 
task  is  equivalent  to  determining  a suitably  weighted  median  of  the  breakpoints. 

In  the  case  of  the  L2-norm,  furthermore,  it  can  be  shown  that  the  vertex-based  iterative 
procedure  converges  to  the  true  optimum  (e.g.  [2]).  That,  however,  does  not  hold  in  the  case 
of  the  Li-norm.  For  this  reason,  two  additional  algorithms  for  Li-optimization  have  been 
developed  in  order  to  improve  the  degree  of  optimality  that  can  be  achieved. 

The  first  of  those  algorithms  is 

edge  — based  elevation  optimization , 

where  the  basic  step  is  to  L\ -optimize  the  elevations  at  the  two  end- vertices  of  a triangulation 
edge  simultaneously  while  keeping  the  elevations  at  all  other  vertices  of  the  triangulation  fixed. 
That  local  optimization  problem  turns  out  to  be  a special  linear  programming  problem,  and 
a special  algorithm  has  been  designed  to  solve  it. 

The  second  algorithm  utilizes  the  fact  that,  for  p G R, 


rv 

(Zp 

~ Zp) 

— 1 - 

1 ' 

II 

1 

— 

K\ 

lrpl 

1 Zp 

Zp  \ 

to  formulate  a 

residual  — weighted  L2  — procedure 

for  Li-approximation.  Roughly  speaking,  we  start  the  ^-approximation  with  equal  weights, 
and  adjust  these  weights  as  we  proceed  on  the  basis  of  observed  residuals.  The  challenge  here 
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is  to  handle  the  cases  of  very  small  residuals  which  result  in  very  large  weights.  There  is  also 
the  question  how  far  to  iterate  the  L2  procedure  before  updating  the  weights.  Our  current 
approach  is  to  adjust  weights  after  one  L2- pass  through  the  vertices  of  the  triangulation. 

Both  algorithms  are  still  suboptimal  and,  as  the  numerical  results  presented  in  this  report 
show,  best  results  so  far  have  been  achieved  by  alternating  both  algorithms.  We  know  that 
these  results  are  suboptimal,  because  it  is  well  known  that,  for  any  fully  optimal  solution  of 
any  Lx  approximation  problem,  the  number  of 

sharp 

elevation  location  points  (fp,  yp)  € R , that  is,  points  for  which  zp  — zp , and  which  are  thus 
represented  exactly  by  the  TIN  surface,  essentially  equals  the  number  of  degrees  of  freedom 
provided  by  the  optimization  parameters.  Here  that  number  is  the  number  |5'\S"|  of  nonfixed 
critical  points.  For  our  best  results,  the  number  of  sharp  points  does  not  quite  reach  that 
goal. 

The  algorithms  and  routines  described  here  were  devised  and  developed  in  the  context  of 
particular  applications  to  terrain  modeling.  For  these  applications,  the  bulk  of  the  elevation 
data  points  were  provided  as  regular  grid  points  with  square  unit  cells  covering  a rectangular 
map  exactly.  A set  S'  of  optional  additional  elevation  points  - “ feature  points ” - were  also 
specified  with  the  proviso  that  they  be  included  into  the  set  of  critical  points  S and  that  their 
given  elevations  are  not  to  be  changed.  If  such  a feature  point  happens  to  be  located  on  the 
grid,  then  its  given  elevation  takes  precedence  over  the  specified  grid  elevation.  Furthermore, 
all  critical  points  other  than  feature  points  had  to  be  grid  points.  This  ensured  that  all 
corresponding  stars  were  complete,  as  defined  above,  because  their  respective  center  vertices 
are  elevation  points.  In  addition,  the  specified  grid  elevations  at  such  critical  points  provided 
good  initial  values  for  iterations. 

For  the  purpose  of  identifying  terrain  artifacts,  it  is  useful  to  generalize  Lx  approximation 
by  scaling  residuals  differently  depending  on  their  sign.  Our  software  thus  includes  an  option 
to  assign  scale  factors  u > 0 if  rp  > 0 and  v > 0 if  rp  < 0.  We  call  that  option 

sign  — oriented  scaling. 

If,  for  instance,  the  object  is  to  identify  artifacts  above  bare  earth,  then  one  would  scale 
positive  residuals  lower  than  negative  residuals.  That  would  reduce  the  influence  of  artifacts 
above  bare  earth  on  the  construction  of  the  TIN  surface.  Those  artifacts  will  then  show  up 
more  clearly  if  the  data  are  compared  to  the  latter.  Sign-oriented  scaling  can  be  specified  for 
any  norm. 

The  algorithms  reported  here  are  as  follows.  For  quantities  not  defined  earlier,  we  ask  the 
reader  to  refer  to  the  detailed  descriptions  provided  later. 

L2FIT(i?,  5,  S',  T,  F,  Z,  Z,  n,  tol) 
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is  a vertex-based  iterative  procedure  for  L2-approximation.  It  is  known  to  actually  converge 
to  the  optimum.  It  has  been  in  use  as  a post-optimization  option  in  the  prototype  NIST  TIN 
package  for  several  years. 

L1FIT1(jR,  S,  S',  T,  F,  Z , Z,  n,  tol , npmax , u,  v ) 

is  a vertex-based  iterative  procedure  for  Lr -approximation  featuring  sign-oriented  scaling.  It 
is  suboptimal.  It  was  the  first  L\  procedure  to  be  considered. 

L1FIT2(.R,  S,  S',  T , F,  E , Z,  Z , n,  tol , npmax , w,  f ) 

is  an  edge-based  iterative  procedure  for  Li-approximation  featuring  sign-oriented  scaling.  It 
is  suboptimal  but  achieves  results  that  are  closer  to  optimality  than  the  vertex-based  version. 
It  was  the  second  L\  procedure  to  be  considered. 

L1FIT3(.R,  S,  S',  T , F , Z,  Z,  n,  tol,  npmax , u,  v,  wscl) 

aims  to  achieve  an  Li  optimum  by  iterating  a residual-weighted  vertex-based  L2  procedure. 
Sign-oriented  scaling  is  included.  The  procedure  is  still  suboptimal  but  by  itself  achieves  a 
higher  degree  of  optimality  than  the  previous  procedure. 

Comparative  numerical  results  are  reported  for  Li-approximation  of  a representative  set 
of  terrain  elevation  data  without  feature  points  and  without  sign-oriented  scaling.  As  was 
mentioned  earlier,  best  results  were  found  when  procedures  L1FIT2  and  L1FIT3  were  alter- 
nated. We  feel  that  those  results  are  close  to  optimal,  particularly,  because  the  number  19558 
of  sharp  points  approaches  the  number  20000  of  critical  points. 


2.  Numerical  results 

The  procedures  presented  here  for  computing  approximate  L\  solutions  have  been  imple- 
mented at  NIST.  They  were  applied  to  a terrain  data  set  consisting  of  a regular  grid  of 
229,761  (521  x 441)  points  with  square  25  x 25 m unit  cells  covering  a rectangular  11  x 13 km 
map.  A triangulation  of  20,000  critical  points  - selected  from  among  the  grid  points  and 
covering  the  map  - was  also  given.  None  of  the  critical  points  were  feature  points  with  fixed 
elevations,  so  S'  = 0.  Sign-oriented  scaling  was  not  used,  so  u — v = 1. 

Table  2.1,  Table  2.2,  and  Table  2.3  show  some  numerical  results  obtained  when  the  vertex- 
based  procedure,  the  edge-based  procedure,  and  the  L2  with  weights  procedure,  respectively, 
were  used  on  the  test  problem.  Table  2.4  shows  some  numerical  results  obtained  when  the 
edge-based  procedure  was  used  on  the  test  problem  after  300  passes  through  the  vertices  by 
the  L2  with  weights  procedure.  In  all  cases,  where  it  applies,  variable  tol  was  set  to  6.058e-04, 
and  variable  wscl  to  1.05.  Given  an  integer  i,  1 < i < n,  where  it  applies,  E(i)  was  set  to  the 
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pass 

CP  Usee 

objective  residual  tolerance 

MAX  residual 

# sharp  points 

1 

121.9 

2508140.342 

158.750 

14976 

2 

121.6 

2466989.055 

158.750 

15015 

5 

121.2 

2461191.750 

155.750 

17990 

10 

121.2 

2460970.594 

155.750 

18604 

20 

121.1 

2460872.057 

155.750 

18638 

30 

121.2 

2460863.037 

155.750 

18639 

50 

121.2 

2460857.512 

155.750 

18637 

75 

121.9 

2460852.957 

155.750 

18639 

100 

121.8 

2460850.886 

155.750 

18640 

125 

121.8 

2460848.991 

155.750 

18644 

150 

121.8 

2460847.192 

155.750 

18649 

Table  2.1:  Results  for  vertex-based  procedure  L1FIT1 


point  in  S which  was  the  first  point  in  the  data  structure  of  the  implementation  that  was  an 
endpoint  of  an  edge  in  T for  which  F(i)  was  the  other  endpoint. 

In  all  tables,  the  first  column  gives  the  sequence  number  of  the  particular  pass.  The  second 
column  displays  its  CPU  running  time  in  seconds.  The  third  column  contains  the  value  of 
the  objective  function  at  the  completion  of  the  pass.  Where  applicable  (Table  2.3),  the  fourth 
column  contains  the  value  of  the  residual  tolerance  used  during  the  pass.  The  fifth  column 
contains  the  maximum  of  the  absolute  values  of  the  residuals  at  the  completion  of  the  pass. 
The  last  column  contains  the  number  of  sharp  grid  points  at  the  completion  of  the  pass,  i.e. 
the  number  of  residuals  whose  absolute  values  were  less  than  0.01  at  the  completion  of  the 
pass. 

We  note  that  the  numerical  results  in  Table  2.1,  Table  2.2,  and  Table  2.3  seem  to  indicate 
that  the  edge-based  procedure  works  better  than  the  vertex-based  procedure,  and  that  the  Z/2 
with  weights  procedure  works  better  than  the  edge-based  procedure.  However,  the  results  in 
Table  2.4  seem  to  indicate  that  using  the  edge-based  procedure  after  the  execution  of  the  L2 
with  weights  procedure  is  a better  approach. 
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pass 

CP  Usee 

objective  residual  tolerance 

MAX  residual 

# sharp  points 

1 

558.4 

2476235.062 

158.750 

15552 

2 

566.4 

2460725.514 

155.750 

16867 

5 

555.0 

2457409.175 

155.750 

18986 

10 

554.6 

2457036.057 

155.750 

19573 

20 

558.8 

2457001.234 

155.750 

19683 

30 

554.9 

2456996.365 

155.750 

19649 

50 

561.0 

2456984.264 

155.750 

19595 

75 

554.5 

2456977.839 

155.750 

19552 

100 

551.4 

2456973.258 

155.750 

19536 

125 

558.9 

2456969.810 

155.750 

19505 

150 

552.5 

2456966.380 

155.750 

19494 

Table  2.2:  Results  for  edge-based  procedure  L1FIT2 


pass 

CP  Usee 

objective 

residual  tolerance 

MAX  residual 

# sharp  points 

1 

121.5 

2588582.203 

146.825 

121.619 

268 

5 

121.5 

2557018.183 

120.793 

126.196 

280 

10 

122.3 

2556963.799 

94.645 

129.640 

278 

20 

122.1 

2554668.834 

58.104 

135.601 

282 

30 

122.1 

2542027.536 

35.671 

137.130 

279 

50 

122.0 

2499697.025 

13.444 

142.532 

262 

75 

122.0 

2469877.300 

3.970 

150.890 

321 

100 

121.9 

2460298.750 

1.172 

154.523 

446 

125 

121.9 

2457317.776 

0.346 

155.564 

900 

150 

121.9 

2456400.229 

1.022e-01 

155.702 

2117 

175 

121.9 

2456122.316 

3.019e-02 

155.737 

6701 

200 

121.9 

2456038.092 

8.915e-03 

155.745 

18046 

225 

121.9 

2456015.632 

2.633e-03 

155.749 

18589 

250 

121.9 

2456013.435 

7.774e-04 

155.749 

18631 

275 

123.1 

2456013.310 

6.058e-04 

155.749 

18603 

300 

122.2 

2456013.128 

6.058e-04 

155.749 

18584 

Table  2.3:  Results  for  L2-based  procedure  L1FIT3 
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pass 

CPUsec 

objective  residual  tolerance 

MAX  residual 

# sharp  points 

1 

545.8 

2455998.827 

155.748 

19161 

2 

553.8 

2455994.372 

155.748 

19318 

5 

537.3 

2455988.608 

155.750 

19545 

10 

539.4 

2455987.781 

155.750 

19580 

20 

539.8 

2455987.481 

155.750 

19576 

30 

539.9 

2455987.377 

155.750 

19578 

50 

547.7 

2455987.331 

155.750 

19570 

70 

540.8 

2455987.276 

155.750 

19569 

150 

539.7 

2455987.265 

155.750 

19558 

Table  2.4:  Results  for  edge-based  L1FIT2  after  L1FIT3 


3.  The  2-dimensional  L\  fitting  problem 

Let  R be  a regular  grid  with  rectangular  cell  units  covering  a rectangular  map  in  the  plane. 
Let  S be  a finite  set  of  critical  points  in  the  map  that  contains  the  four  corners  of  the  map, 
and  let  T be  a triangulation  for  S.  Let  S'  be  the  set  of  feature  points,  that  is,  of  critical 
points  with  prescribed  fixed  elevations.  For  each  grid  point  p in  R assume  that  an  elevation 
zp  is  associated  with  p.  Given  p in  R , let  f be  a triangle  in  T that  contains  p,  and  define 
Si(p),  i = 1,2,3,  as  the  points  in  S that  are  the  vertices  of  t.  Given  p in  R , define  A i(p), 
i = 1,2,3,  as  the  barycentric  coordinates  of  p relative  to  T so  that  p equals  \{p)si{p) 
and  0 < Xl(p)  < 1 for  each  i,  i = 1, 2, 3.  For  each  point  s in  S'  assume  that  an  elevation  zs  is 
associated  with  s.  We  search  then  for  a set  of  elevations  {zs  : s £ S \ S'}  so  that  the  following 
function  is  minimized: 

3 

/G  I Zp  ~ ^2  ^i^P)ZSi(p)  I- 

p£/?  i—  1 

More  generally,  given  u , v > 0,  and  if  for  each  p in  R we  let  zp  equal  Y2=i  ^ i(p)zsl(p p we 
search  then  for  a set  of  elevations  {zs  : s G S \ 5'}  so  that  the  following  function  is  minimized: 

(3.i)  « E E 

PG.R+  peR~ 

where  R+  = {p  G R : zp  — zp  > 0}  and  R~  = {p  E R : zp  — zp  < 0}. 
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4.  The  2-dimensional  L2  fitting  problem 

Of  related  interest  is  the  quadratic  version  of  the  problem.  Here  we  search  for  a set  of  elevations 
{:s  : s € 5 \ 5'}  so  that  the  following  function  is  minimized: 

3 

(4.1)  ^2  (ZP  — ^2  ^i(p)zsi{p))  ■ 

p£R  i=  1 

Since  this  function  is  quadratic  it  has  a unique  minimum  point  which  is  the  point  at  which 
its  gradient  equals  zero  [2]. 

Given  a point  s in  S we  define  the  star  of  s as  the  union  of  the  triangles  in  T that  have  s 
as  a vertex,  and  denote  by  Rs  the  set  of  grid  points  that  it  contains. 

For  a fixed  point  s in  S’,  we  assume  without  any  loss  of  generality  that  for  each  point  p in 
Rs,  si(p)  equals  s.  Under  this  assumption  the  partial  derivative  of  Function  4.1  with  respect 

to  zs  is 

2 (ZP  ~ A1  (p)zs  - A2 {p)zs2(p)  - A3 {p)zS3(p)){-^l{p))- 

pERs 

Setting  this  partial  derivative  to  zero  we  are  then  able  to  solve  for  2S  and  obtain 

z8  = {'52{zp~  A2  {p)zS2(P)  ~ h{p)zs3(P))^i{p))/  J2  ( Xi(P ))2- 

p£Rs  pERs 

The  following  iterative  procedure,  which  is  based  on  the  above  formula,  can  be  used  for 
computing  the  minimum  point  of  Function  4.1.  That  it  terminates  in  a finite  number  of  steps 
follows  again  from  the  fact  that  Function  4.1  is  quadratic  [2],  Here  we  assume  that  there  are 
n points  in  5,  n a positive  integer,  and  that  the  n points  are  ordered  in  some  fashion.  A 
one-to-one  function  F from  {l,...,n}  onto  S is  defined  by  setting  F(i)  to  the  ith  point  in 
S for  each  i,  i = 1, . . . , n.  A real- valued  elevation  function  Z with  domain  S is  defined  by 
setting  Z(s)  to  2S  for  each  s in  S',  and  to  zs  otherwise.  At  the  end  of  the  execution  of  the 
procedure  Z will  contain  the  minimum  point.  Another  real- valued  elevation  function  Z with 
domain  R is  defined  by  setting  Z(p)  to  zp  for  each  p in  R.  The  procedure,  called  L2FIT, 
requires  an  input  variable  called  tol  which  must  be  set  to  a positive  number  and  is  used  as 
a tolerance  or  adjustment  constant  during  the  execution  of  the  procedure.  Two  procedures, 
called  STARGRIDPOINTS  and  BARYCENTRIC  are  used  as  primitives  in  L2FIT.  Given 
integer  i,  1 < i < n,  with  F(i)  G S \ 5',  for  some  positive  integer  m,  STARGRIDPOINTS 
locates  points  G(j),  j = 1, . . . , m,  that  are  the  points  in  Rp(i).  Given  an  integer  j,  1 < j < m, 
BARYCENTRIC  locates  a triangle  in  T with  F(i)  as  a vertex  that  contains  G(j),  identifies 
the  other  two  vertices  Q2,  Q 3 of  the  triangle,  and  computes  the  barycentric  coordinates  Rfi, 
W2,  IU3  of  G(j)  relative  to  the  triangle  so  that  G(j)  equals  W\  ■ F(i)  + W2  ■ Q2  + W3  ■ Q3. 
The  outline  of  the  procedure  follows. 
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procedure  L2FIT(R,  S,  S' , T,  F,  Z,  Z,  n,  tol ) 
begin 
flag  :=  0; 

while  (flag  = 0)  do 
begin 
flag  :=  1; 

for  z :=  1 until  n do 
begin 

if  (F(z)  E S \ S')  then 
begin 

(G,  m):=STARGRIDPOINTS(F(«),  R,  T); 
numerator  :=  0.0;  denominator  :=  0.0; 

for  j :=  1 until  m do 
begin 

(Q2,  Qs,  WY,  VF2,  VF3):=BARYCENTRIC(T,  G(j),  F(t)); 
numerator  :=  numerator + 

(Z(G(j))  - W2  ■ Z(Q2)  - W3  ■ Z(Qs))  ■ Wu 
denominator  :=  denominator  + (VFi)2 

end 

zneie  :=  numerator / denominator] 
if  (I  Z(F(i))  — znew  |>  tol ) then 
begin 

flag  :=  0; 

Z(F(i))  :=  zneu> 

endendendendend 


5.  Vertex-based  iterative  procedure  for  L\  approximation 

An  approach  similar  to  the  one  in  the  previous  section  can  be  used  for  obtaining  approximate 
solutions  to  L\  fitting  problems.  Given  s in  S \ Sf  for  each  point  s in  S \ {s}  we  assume  that 
an  elevation  is  associated  with  s.  We  search  then  for  an  elevation  so  that  the  function 

3 

yi  i % ~ \ (p)zsi(p ) i 

pGR  i= 1 


is  minimized. 
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More  generally,  given  u , v > 0,  we  search  for  an  elevation  so  that  the  function 

3 3 

U Z (ZP  ~ Z Xl{P)ZSi(p))  -V  (M  ~ E Xi(P)ZSr(p )) 
p<Etf+  J=1  pGR-  i=  1 


is  minimized,  where 

i?+  = {p  € R : ip  - Ei=i  Xi{p)zsi(P)  > 0}  and  R~  = {p  G R : zp  - £-= i Ai(p)25i(p)  < 0}. 

Assume  without  any  loss  of  generality  that  for  each  point  79  in  Rs,  Si(p)  equals  s,  and  for 
some  positive  integer  m,  let  pj,  j = 1, . . . , m be  the  points  in  Rs.  The  above  problem  is  then 
equivalent  to  searching  for  an  elevation  zs  so  that  the  function 

m 

(5.1)  Y.  Wj\Zpj  - A1  (Pj)zs  - A2 {Pj)z82{pj)  - X3(Pj)zS3(Pj)\ 

j- 1 

is  minimized,  where  Wj  equals  u if  zPj  — Xi(pj)zs  — A 2(Pj)zS2(Pj)  — A3 (Pj)zS3(Pl)  > 0,  and  v oth- 
erwise, for  each  j,  j = 1, . . . ,m.  This  function  is  convex,  continuous,  piece-wise  linear,  and 
has  at  least  one  minimum  point. 

We  show  how  a minimum  point  of  Function  5.1  can  be  found.  For  each  j , j = 1, . . . ,m, 
assume  without  any  loss  of  generality  that  Xi(pj)  is  not  equal  to  zero,  and  define  c3  as  follows: 

Cj  = (zPj  - X 2{Pj)zS2{pj)  - A3 (py ) ) ) / A 1 (Py ) - 

Assume  without  any  loss  of  generality  that  m > 1,  and  that  for  each  j,  j = 1, . . . , m — 1, 
Cj  < Cj+i . For  each  j,  j = 1, . . . , m — 1,  let  Ij  be  the  interval  [cj,Cj+i\.  Define  /0  as  the 
interval  (— 00,  Ci],  and  Im  as  the  interval  [cm,oo).  In  the  interval  /0  the  slope  of  Function  5.1 
is  — u EjfcLi  Xi(pk),  and  in  Im  it  is  v E£Li  Ai (Pfe ) - For  each  j,  j = 1, . . . , m — 1,  in  the  interval 
Ij  the  slope  is  v Ei=i  Ai(Pfc)  — u Y1T=j+i  Xi{pk)-  Clearly,  the  slope  must  change  signs  at  an 
endpoint  of  one  the  intervals,  and  this  is  then  a minimum  point  of  Function  5.1. 

The  following  iterative  procedure  for  attempting  to  minimize  approximately  Function  3.1, 
is  based  on  the  ideas  presented  above  for  minimizing  Function  5.1.  The  procedure,  called 
L1FIT1,  requires  an  input  variable  called  npmax  which  must  be  set  to  a positive  integer 
and  is  used  as  the  maximum  number  of  passes  allowed  through  the  vertices,  thus  guaran- 
teeing that  the  execution  of  the  procedure  terminates.  Besides  STARGRIDPOINTS  and 
BARYCENTRIC,  a procedure  called  INCREASORT  is  used  as  a primitive  in  L1FIT1.  Given 
an  integer  mw , mw  > 1,  a one-to-one  function  M from  {1, . . . ,mw}  onto  a subset  / of  the 
positive  integers,  and  a function  C from  I into  the  set  of  real  numbers,  INCREASORT  rear- 
ranges M so  that  for  each  j,  j — 1, . . . , mw  — 1,  C(M(j ))  < C(M ( j + 1)).  The  outline  of  the 
procedure  follows. 
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procedure  LlFITl(i?,  S , S',  T,  F,  Z , Z,  n,  fo/,  npmax,  u , n) 

begin 

flag  0;  npass  :=  0;  w :=  u + v, 
while  {f  lag  = 0 and  npass  < npmax ) do 

begin 

flag  1;  npass  npass  + 1; 

for  i 1 until  n do 
begin 

if  {F{i)  e S \ S')  then 
begin 

(G,  ra):=STARGRIDPOINTS(F(z),  R,  T); 

mw  :=  0;  slope  :=  0.0; 

for  j :=  1 until  m do 
begin 

(Qa,  Q3,  Wi,  W2,  IE3):=BARYCENTRIC(T,  G(j),  F(z)); 

if  (MR  > 0.0)  then 
begin 

mw  :=  mw  + 1; 

C(mW)  :=  (Z(G(i))  - W2  ■ Z(Q2 ) - 1P3  • ^(Os))/^ 
D(mw ) :=  MA; 

M {m,w)  mw ; 

slope  :=  slope  — u ■ W\ 

endend 

if  (raw  > 1)  then  INCREASORT(M,  C,  raw); 

J :=  0; 

while  ( slope  < 0.0)  do 

begin 

j :=  j + 1; 

slope  :=  s/ope  + w • D(M(j )) 

end 

znew  :=  C(M(j )); 
if(|Z(F(t))  — 2:new  |>  tol ) then 
begin 
flag  :=  0; 

Z(F(i))  znew 
endendendendend 
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6.  Edge-based  iterative  procedure  for  L\  approximation 

Given  s,  s'  in  S \ S',  s and  s'  the  endpoints  of  an  edge  of  a triangle  in  T,  and  assuming  that 
an  elevation  z$  is  associated  with  each  point  s in  S \ {.s,  s'},  we  now  search  for  elevations 
and  zs>  so  that  the  function 

3 3 

u {Zp  ~ yi  \ {p)zsi(v))  — v {zp  — ^j{p) zsj(p)) 

peR+  i= 1 peR - i= 1 

is  minimized.  Here  R+,  R~,  u,  v are  as  in  the  previous  section. 

Assume  without  any  loss  of  generality  that  for  each  point  p in  Rs , Si(p)  equals  s,  that  for 
each  point  p in  Rs>,  S2{p)  equals  s',  and  that  the  sets  Rs  \RS>  and  Rs'  \RS  are  nonempty. 
For  positive  integers  mi,  m2,  m3,  let  pji,  j = be  the  points  in  RsnRs',  let  pj2, 

j = 1, . . . , m2  be  the  points  in  Rs  \ Rs',  and  let  p3 3,  j — 1, . . . , m3  be  the  points  in  Rs>  \ Rs. 
The  above  problem  is  then  equivalent  to  searching  for  elevations  zs  and  zs>  so  that  the  following 
function  is  minimized: 

mi 

(6-1)  J2wij\2Pn  ~ xi(Pn)zs  ~ X2(P]i)zs'  ~ h(Pji)zs3(p3l)\  + 

3= 1 
m2 

'52w2j\zPj2  - Ai  {Pj2)zs  - X2(Pj2)zS2(pl2)  - A3(pj2)^S3(pi2)|  + 

3=1 

m3 

w3j  I zPj3  — Ai  (Pj3)2Si(pj3)  — A2(pj3)2s/  — A3(pj3)2;S3(pj3)  |. 

Here  wkj  equals  u if  zPjk  - xi{Pjk)zSl(Pjk)  - A 2(Pjk)zS2(Pjk)  - xs{Pjk)zS3{Pjk)  > 0,  and  v other- 
wise, for  each  k,  j,  k = 1, 2,  3,  j — 1, . . . , m*,.  In  what  follows  we  think  of  Function  6.1  as  a two- 
variable  real- valued  function  whose  first  variable  is  zs>  and  whose  second  variable  is  zs.  We  de- 
note it  by  / so  that  if  and  zs>  are  set  to  values,  f(zs> , zs)  then  denotes  the  corresponding  value 
of  the  function.  Finally,  we  denote  by  H(s ',  s)  the  set  of  triplets  of  the  form  (zs>,  zs,  f(zs>,  zs)). 
H(s',  s)  is  a convex,  continuous,  piece-wise  linear  surface  in  3-dimensional  space,  and  zs,  zs> 
exist  at  which  / achieves  its  minimum  value,  and  (zs>,  zs,  f{zsi,  zs))  is  a vertex  of  H(sr,  s). 

For  each  j,  j = 1, . . . , mi,  assume  without  any  loss  of  generality  that  Ai(pji)  is  not  equal 
to  zero,  and  define  Cj,  bj , d3  as  follows: 

cj  — [zp.  1 — A3(pji)zS3(Pjl))/ Ai(pji), 

bj  = -A2(pji)/A1(pjl), 
dj  = X\{pji). 
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For  each  j,  j = 1 , . . . , m2,  assume  without  any  loss  of  generality  that  A i(pj2)  is  not  equal 
to  zero,  and  define  cJ+Tni,  h,+mi,  dj+mi  as  follows: 

Cj+mi  = ( zPj2  ~ ^(Pj2)zs2{pj2)  ~ ^3{‘Pj2)zS3(pj2)) / (Pj?)’ 
bj+rrii  ~ O-Oj 
dj+rrii  = ^l\Pj2) • 

For  each  j,  j = 1, . . . , m3,  assume  without  any  loss  of  generality  that  A2(pj3)  is  not  equal 
to  zero,  and  define  c' , d ' as  follows: 

Cj  = (%j3  — ^1  (Pj3)zsi(pj3)  — A 3 (pj  3 ) Zg3  (pj3  ) ) / A2  (Pj3  ) , 
d'j  = A2(pj3). 

Setting  mru  to  m^  + m2,  and  mv  to  m3,  for  each  j,  j = denote  by  lj  the 

straight  line  in  the  2V  — plane  defined  by  the  linear  equation  zs  — bjzsi  + Cj,  and  for  each  j, 
j = 1, . . . , mu,  denote  by  /'  the  straight  line  in  the  same  plane  defined  by  the  linear  equation 
z's  = dj . Clearly,  from  the  polyhedral  nature  of  H(s',s)  it  follows  that  the  perpendicular 
projections  onto  the  zs<  — zs  plane  of  the  1-dimensional  faces  of  H(s s)  are  contained  in  these 
lines,  and  those  of  the  vertices  are  points  where  these  lines  intersect.  Thus,  the  problem  of 
minimizing  Function  6.1  reduces  to  that  of  finding  among  the  points  at  which  at  least  two  of 
the  lines  lj , j — 1, . . . , mw , , j = 1, . . . , mu,  intersect,  one  at  which  the  function  achieves  its 

smallest  value. 

Given  an  integer  mopt,  1 < mopt  < mw , we  show  how  a minimum  point  of  /,  i.  e. 
Function  6.1,  restricted  to  line  lmopt  can  be  found  in  a manner  similar  to  the  one  used  in  the 
previous  section  for  minimizing  Function  5.1.  Without  any  loss  of  generality  we  assume  that 
at  any  point  in  the  zs<  — zs  plane  at  most  two  of  the  lines  lJ7  j = 1, . . . , mw,  l'rj,  j = 1, . . . , mu, 
intersect.  Thus,  if  two  of  these  lines  intersect,  and  Function  6.1  restricted  to  the  two  lines 
achieves  its  smallest  value  at  the  point  at  which  the  two  lines  intersect,  then  from  the  convexity 
and  the  polyhedral  nature  of  H(s',s ) it  follows  that  among  the  points  at  which  any  two  of 
the  lines  intersect,  at  this  point  the  function  achieves  its  smallest  value.  Therefore,  as  pointed 
above,  this  then  implies  that  at  this  point  the  function  also  achieves  its  minimum  value  without 
any  restrictions  on  zs  and  zs> . 

For  each  j,  j = 1, . . . ,mw , j ^ mopt , assume  without  any  loss  of  generality  that  l0  is  not 
parallel  to  lmopti  and  define  a3  as  follows: 


Uj  — ( Cj  (-mopt ) / ( bmopt  bj  ) . 


Clearly,  aj  is  the  value  of  zst  at  which  lj  intersects  lm0pt ■ Assume  without  any  loss  of 
generality  that  mw  > 2,  that  mopt  equals  mw,  and  that  for  each  j , j = l,...,mw  — 2, 


15 


dj  < aJ+ 1-  For  the  sake  of  simplicity  and  only  for  the  current  discussion,  assume  that  mv 
equals  zero,  and  for  each  j,  j = 1, . . . ,mw  — 2,  let  Ij  be  the  interval  [a,j,aj+ 1]  in  the  zs' 
axis  of  the  zs‘  — zs  plane.  Define  /0  as  the  interval  (— 00,0^,  and  Imw-\  as  the  interval 
[«„IU,_i,oo)  in  the  same  axis.  Define  a one-variable  real-valued  function  f by  setting  f'{zs>) 
to  f {zs> , bmoptzst  + cmopt)  for  each  value  of  zs>.  f is  convex,  continuous,  piece- wise  linear, 
and  has  at  least  one  minimum  point.  For  each  j,  set  Uj  to  u if  bmopt  — bj  is  positive,  and 
to  — v otherwise,  and  set  Vj  to  v if  bmopt  — bj  is  positive,  and  to  — u otherwise.  With  this 
terminology,  in  the  interval  I0  the  slope  of  f is  — uj{bm0pt  — bj)dj,  and  in  Imw-i 

it  is  T!k=\l  vj(bmopt  ~ bj)dj.  For  each  j , j = 1, . . . ,mw  - 2,  in  the  interval  Ij  the  slope  is 
aliMw-w-as;!  Uj  {bmopt  ~ bj)dj.  Clearly,  the  slope  of  /'  must  change  signs  at 
an  endpoint  of  one  of  the  intervals,  and  this  is  then  a minimum  point  of  Thus,  given  zs> 
at  which  /'  achieves  its  minimum  value  it  follows  that  Function  6.1  restricted  to  lm0pt  has  a 
minimum  point  at  {zs>,  bmoptzs’  + cmopt ). 

A minimum  point  of  Function  6.1  without  any  restrictions  on  zs  and  zs<  can  then  be  found 
through  a two-step  iterative  procedure  as  follows: 

Step  1.  Mark  each  one  of  the  lines  Ij , j = 1 Ij,  j = 1 ,...,mv  as  not  optimized. 

Select  arbitrarily  one  of  these  lines  and  call  it  the  current  line. 

Step  2.  Find  a minimum  point  of  Function  6.1  restricted  to  the  current  line,  and  call  this 
point  the  current  minimum  point.  Mark  the  current  line  as  optimized , identify  the  one 
line  among  the  lines  Ij,  j = 1, , mw , , j — 1, ... , mv,  that  intersects  the  current 
line  at  the  current  minimum  point,  and  call  this  line  the  next  line.  If  the  next  line 
has  been  marked  as  optimized  then  the  current  minimum  point  is  a minimum  point  of 
Function  6.1  without  restrictions  and  the  procedure  terminates.  Else  call  the  next  line 
the  current  line  and  go  back  to  the  beginning  of  Step  2. 

The  following  iterative  procedure  for  attempting  to  minimize  approximately  Function  3.1, 
is  based  on  the  ideas  presented  above  for  minimizing  Function  6.1.  Here  we  assume  that 
for  each  i,  i — 1, . . . , n,  an  edge  of  a triangle  in  T is  selected  so  that  F{i)  is  an  endpoint  of 
this  edge.  A function  E from  {1, . . . , n}  into  S is  defined  by  setting  for  each  i , i — 1, . . . , n, 
E{i)  to  the  point  in  S which  is  the  other  endpoint  of  the  edge  selected  with  F{i ) as  an 
endpoint.  Besides  INCREASORT,  the  procedure,  called  L1FIT2,  requires  primitives  STAR- 
GRIDPOINTS1,  STARGRIDPOINTS2,  STARGRIDPOINTS3  as  well  as  BARYCENTRIC1, 
B ARY  CENTRIC2,  and  BARYCENTRIC3.  Given  integer  i,  1 < i < n,  with  F{i)  eS\S',  for 
some  positive  integer  mi,  STARGRIDPOINTS1  locates  points  G\{j),  j = 1, . . . , m1?  that  are 
the  points  in  Rf^)  fl  Re{i)-  Given  that  Rf(i)  \ Re{i)  is  not  empty,  for  some  positive  integer  m2, 
STARGRIDPOINTS2  locates  points  G2{j ),  j — 1,  • • • , nr,2,  that  are  the  points  in  Re(i)  \ Re{i)- 
It  sets  m2  to  zero  if  i?F(i)  \ Re(i)  is  empty.  Given  that  Re(i)  \ Rf{i)  is  not  empty,  for  some  posi- 
tive integer  m3,  STARGRIDPOINTS3  locates  points  G3(j),  j = 1, . . . , m3,  that  are  the  points 
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in  Re(i)  \ -Rf(i)-  It  sets  m3  to  zero  if  Re(i)  \ Rf{i)  is  empty.  Given  an  integer  j,  1 < j < m i, 
BARYCENTRIC1  locates  a triangle  in  T with  F(i)  and  E(i)  as  vertices  that  contains  G\(j), 
identifies  the  third  vertex  Q3  of  the  triangle,  and  computes  the  barycentric  coordinates  Wi, 
W2,  W3  of  G\{j)  relative  to  the  triangle  so  that  G i (j ) equals  W\  ■ F(i)  + W2  • E(i)  + W3  • Q3. 
Assuming  m2  is  positive,  given  an  integer  j,  1 < j < ra2,  BARYCENTRIC2  locates  a triangle 
in  T with  F(i)  as  a vertex  that  contains  G2(j ),  identifies  the  other  two  vertices  Q2 , Q3  of 
the  triangle,  and  computes  the  barycentric  coordinates  Wi,  W2,  W3  of  G2(j)  relative  to  the 
triangle  so  that  G2 (j ) equals  W\  • F(i)  + W2  ■ Q2  4-  W3  ■ Q3.  Finally,  assuming  m3  is  positive, 
given  an  integer  j,  1 < j < m3,  BARYCENTRIC3  locates  a triangle  in  T with  E(i)  as  a vertex 
that  contains  G3 (j ) , identifies  the  other  two  vertices  Qi,  Q3  of  the  triangle,  and  computes 
the  barycentric  coordinates  W3,  W2,  W3  of  G3(j)  relative  to  the  triangle  so  that  G3  (j ) equals 
H i ' Q\  + ■ E(i)  + W j • Q3. 

The  outline  of  the  procedure  follows.  In  lines  6-41  of  the  procedure,  the  nonvertical 
straight  lines  (lj,  j = l,...,raie,  above),  and  the  vertical  straight  lines  (/',  j — 1 
above)  associated  with  the  current  edge  are  identified.  In  lines  42-60  of  the  procedure  the 
first  straight  line  to  be  called  the  current  line  (for  the  current  edge)  is  identified,  and  it  is  the 
one  nonvertical  line  at  which  Function  6.1  (for  the  current  edge)  attains  its  minimum  for  zs< 
arbitrarily  close  to  — oo.  In  line  66  the  zs>  values  where  the  vertical  lines  intersect  the  zs/-axis 
are  sorted  in  increasing  order.  The  zs>  values  where  the  current  line  intersects  nonvertical  lines 
are  computed  in  lines  72-82  of  the  procedure  if  the  current  line  is  nonvertical,  and  in  lines 
111-114  if  it  is  vertical.  These  values  are  sorted  in  increasing  order  in  line  84  if  the  current  line 
is  nonvertical,  and  in  line  115  if  it  is  vertical.  In  lines  87-99  if  the  current  line  is  nonvertical, 
and  in  lines  117-120  if  it  is  vertical,  a straight  line  to  be  called  the  next  line  is  identified  which 
is  either  a nonvertical  or  a vertical  line  whose  intersection  with  the  current  line  is  a minimum 
point  of  Function  6.1  restricted  to  the  current  line  (to  be  called  the  current  minimum  point). 
The  current  line  is  marked  as  optimized  in  line  106  of  the  procedure  if  it  is  nonvertical  and 
in  line  109  if  it  is  vertical.  Given  that  the  next  line  has  already  been  marked  as  optimized, 
in  lines  122-128  of  the  procedure  if  the  current  line  is  vertical,  and  in  lines  129-134  if  it  is 
nonvertical,  another  nonvertical  line,  if  any,  to  be  called  the  current  line  is  identified  which  is 
not  marked  as  optimized  and  that  contains  the  current  minimum  point.  If  no  such  a line  exists 
then  the  current  minimum  point  is  a minimum  point  for  Function  6.1  without  restrictions  and 
if  necessary  the  values  of  and  zs>  (for  the  current  edge)  are  corrected  in  lines  135-140  of  the 
procedure. 

procedure  L1FIT2(R,  A,  S' , T,  T,  E , Z,  Z,  n,  tol , npmax , u , v ) 

begin 

1.  flag  1 :=  0;  npass  0;  w :=  u + v; 

2.  while  (flag  1 = 0 and  npass  < npmax ) do 

begin 
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3. 

4. 


6. 
i . 
8. 

9. 

10. 
11. 

12. 

13. 

14. 

15. 

16. 

17. 

18. 

19. 

20. 
21. 
22. 

23. 

24. 

25. 

26. 

27. 

28. 
29. 


flag  1 :=  1;  npass  :=  npass  + 1; 

for  i :=  1 until  n do 
begin 

if  (F(i)  G S \ S')  then 
begin 

{Gu  mi ) :=STARGRIDP0INTS1  (F(i) , E(i),  R,  T); 

(G2,  m2):=STARGRIDPOINTS2(F(z),  E(i),  R,  T); 
if  {E{i)  eS\S')  then  (G3,  m3):— STARGRIDPOINTS3(F(?),  E(i),  R,  T) 
else  m3  :=  0; 
mw  :=  0;  mv  :=  0; 
for  j :=  1 until  m\  do 
begin 

(Q3,  Wuw2,  W'3):=BARYCENTRIC1(T,  G,  (j).  F(i),  £(»)); 
if  (RA  > 0.0  and  E'(z)  G S \ S')  then 
begin 

mw  :=  mw  4-  1; 

C(mw)  :=  ■ Z(Q3))/Wp, 

B(mw)  :=  - W2/Wi; 

D(mw)  :=  Wi 

end 

elseif  (Wi  >0.0)  then 
begin 

mw  :=  mw  + 1; 

C(mw)  :=  (Z(Gi(j))  - W2  ■ Z(E(i))  - W3  ■ Z(Q3))/W{ ; 

B(mw)  :=  0.0; 

D(mw)  :=  IPi 

end 

elseif  (W2  > 0.0  and  E(i)  e S \ S')  then 
begin 

mv  :=  mv  + 1; 

C'(mv)  :=  (Z(G1(j))  - IV,  • Z(Q3))/W2 ; 

D’imv)  :=  W2 

endend 

if  (m2  > 0)  then 
begin 

for  j :=  1 until  m2  do 
begin 

(O2,  Os,  WUW2,  W3):=BARYCENTRIC2(T,  G2(i),  F(i)); 
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30. 

31. 

32. 

33. 

34. 

35. 

36. 

37. 

38. 

39. 

40. 

41. 

42. 

43. 

44. 

45. 

46. 

47. 

48. 

49. 

50. 

51. 

52. 

53. 

54. 


if  (Wi  > 0.0)  then 
begin 

mw  :=  mw  + 1; 

C(mw)  :=  (Z(G2(j))  - W2  • Z(Q2)  - W3  ■ Z(Q3))/Wi; 
B(mw)  :=  0.0; 

D(mw)  :=  W\ 

endendend 
if  (m3  > 0)  then 
begin 

for  j :=  1 until  m3  do 
begin 

(Qi,  Qa,  WuW2,  fy3):=BARYCENTRIC3(T,  G3(j), 
if  (W2  > 0.0)  then 
begin 

mv  :=  mv  + 1; 

C’(mv ) :=  (Z(G3(j))  ~ W,  ■ Z(Q ,)  - W3  ■ Z(Q3))/W2; 
D'(mv)  :=  W2 

endendend 

slopew  :=  0.0; 

for  j :=  1 until  mw  do 
begin 

M{j)  ■■=  j\  NU)  :=  0; 

slopew  :=  slopew  + v • D(j) 

end 

if  (mw  > 1)  then 
begin 

INCREASORT (M,  B,  mw ); 

for  j :=  1 until  mw  — 1 do 
begin 

33  ■=  3 + !; 

while  (jj  < mw  and  B(M(jj ))  = B(M(j )))  do 

begin 

if  (C(M(jj))  > C(M(j)))  then 

begin 

mt  :=  M(j ); 

M(j)  :=  M(jj); 

M(jj)  :=  mt 

end 
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55. 

jj  ■=  jj  + 1 

endendend 

5G. 

j ■■=  0; 

57. 

while  ( slopew  > 0.0)  do 

begin 

58. 

j j + 1; 

59. 

slopew  slopew  — w ■ D(M(j)) 

end 

60. 

mopt  M(j ); 

61. 

slopev  :=  0.0; 

62. 

if  ( mv  > 0)  then 

begin 

63. 

for  j :=  1 until  mv  do 

begin 

64. 

o 

II 

■'-s 

'^~s 

II 

■'-J 

65. 

slopev  :=  slopev  + u • D'(j ) 

endend 

66. 

if  (mv  > 1)  then  INCREASORT 

67. 

znwa  :=  C(mopt );  znwb  :=  Z(E(i))\ 

68. 

flagl  :=  0; 

69. 

while  ( flag2  = 0)  do 

begin 

70. 

flag2  :=  1;  mb  0; 

71. 

mbc  :=  1;  m.vc  :=  1;  slopew  :=  slopev 

72. 

for  j :=  1 until  mw  do 

begin 

73. 

b :=  B(mopt)  — B(j ); 

74. 

if  (b  > 0.0)  then 

75. 

slopew  :=  slopew  + u ■ b ■ D(j ) 

76. 

elseif  (b  < 0.0)  then 

77. 

slopew  :=  slopew  — v ■ b • D(j). 

78. 

if  (b  ^ 0.0)  then 

begin 

79. 

c :=  C(j)  - C (mopt) ; 

80. 

mb  :=  mb  + 1; 

81. 

M(mb)  :=  j; 

82. 

^4 O')  :=  c/b 

endend 

83. 

if  (m6  > 0 or  mv  > 0)  then 
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84. 

begin 

if  ( mb  > 1)  then  INCREASORT(M,  A,  mb); 

85. 

M (mb  + 1)  :=  mw  + 1;  A(mw  + 1)  :=  0.0; 

86. 

M'(mv  + 1)  :=  mv  + 1;  C'(mv  + 1)  :=  0.0; 

87. 

while  ( slopew  > 0.0)  do 

88. 

begin 

if  (( mbc  < mb  and  mvc  < mv  and  A(M (mbc)) 

89. 

or  ( mbc  < mb  and  mvc  > mv))  then 

begin 

flagS  :=  1; 

90. 

jb  M (mbc); 

91. 

b :=  B(mopt)  — B(jb); 

92. 

if  (b  > 0.0)  then 

93. 

slopew  slopew  — w • b ■ D(jb ) 

94. 

else 

slopew  :=  slopew  + w • b • D(jb); 

95. 

mbc  :=  mbc  + 1 

96. 

end 

else 

begin 

flag3  :=  2; 

97. 

jv  :=  M'(mvc); 

98. 

slopew  slopew  — w ■ D'(jv); 

99. 

mvc  :=  mvc  + 1 

100. 

endend 

if  ( flag3  = 1)  then 

101. 

begin 

znwb  :=  A(jb); 

102. 

if  (N(jb)  ^ 0)  then  flag3  :=  3 

103. 

end 

else 

begin 

znwb  :=  C'(jv); 

104. 

if  (N'(jv)  ^ 0)  then  flag3  :=  4 

105. 

end 

znwa  :=  B(mopt)  ■ znwb  + C(mopt); 

106. 

N(mopt)  1; 

107. 

if  (flag3  = 1)  then  mopt  :=  jb;  ,flag2  :=  0 

108. 

elseif  (flag3  = 2)  then 
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begin 


109. 

N’(jv)  :=  1; 

110. 

slopew  0.0; 

111. 

for  j 1 until  mw  do 

begin 

112. 

M(j)  :=  j; 

113. 

slopew  :=  slopew  + u • D(j); 

114. 

A(j)  :=  B(j)  • znwb  + C(j) 

end 

115. 

if  ( mw  > 1)  then  INCREASORT(il7,  A,  mw)] 

116. 

J :=  0; 

117. 

while  ( slopew  > 0.0)  do 

begin 

118. 

j ■—  j + 1; 

119. 

slopew  :=  slopew  — w ■ D(M(j )) 

end 

120. 

jb  :=  MU); 

121. 

if  ( N(jb ) = 0)  then  mopt  :=  jb ; flag2  0 

else 

begin 

122. 

zmm  :=  A(jb ); 

123. 

6 :=  0.0;  k :=  1;  jj  :=  j + 1; 

124. 

while  (6  = 0.0  and  k = 1 and  jj  < mw)  do 

begin 

125. 

b :=  A(M(jj ))  — znwa\ 

126. 

k :=  7V(M(jj )); 

127. 

JJ  ■=  JJ  + 1 

end 

128. 

if  (5  = 0.0  and  k = 0)  then  mopt  :=  M (jj  — 

endend 

else 

begin 

129. 

6 :=  0.0;  k :=  1;  jj  :=  m6c; 

130. 

while  (6  = 0.0  and  A;  = 1 and  jj  < mb)  do 

begin 

131. 

b :=  A(M(jj ))  — znwb ; 

132. 

* :=  N(MUj)); 

133. 

jj  :=  jj  + 1 

end 

; flag 2 
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134.  if  ( b = 0.0  and  k = 0)  then  mopt  :=  M(jj  — 1);  flag 2 :=  0 

endendend 

135.  if  (|  Z(F(i))  — znwa  |>  tol)  then 

begin 

136.  flagl  :=  0; 

137.  Z(F(i))  :=  znwa 

end 

138.  if  (|  Z(E(i))  — znwb  |>  tol ) then 

begin 

139.  flagl  :=  0; 

140.  Z(E(i))  :=  zmu& 

endendendendend 


7.  Residual-weighted  L2  iterative  procedure  for  L\  approximation 

An  approach  based  on  the  L2  iterative  procedure  presented  in  Section  4 can  be  used  for 
obtaining  approximate  solutions  to  Li  fitting  problems.  This  is  based  on  the  fact  that  given 
a real  number  r if  r does  not  equal  zero  then  |r|  equals  r2/|r|. 

Given  s in  S\S' , we  assume  that  an  elevation  z and  a positive  number  wcut  are  associated 
with  s.  In  addition,  for  each  point  s in  S \ {s}  we  assume  that  an  elevation  z~s  is  associated 
with  s.  Given  u,  v > 0,  and  assuming  again  without  any  loss  of  generality  that  for  each  point 
p in  Rs,  Si(p)  equals  s,  and  that  for  some  positive  integer  m,  pj,  j = 1, . . . ,m,  are  the  points 
in  Rs,  we  search  then  for  an  elevation  so  that  the  function 

m 

(7-1)  ^2(wj/wj){zp.  - X i{pj)zs  - X2(pj)zS2iPj)  - X3{Pj)zS3(Pj))2 

3= 1 

is  minimized.  Here  Wj  equals  u if  zPj  — X\ (pj)z  — X2 {Pj)zS2(pg  — \3(pj)zS3(p  ) > 0,  and  v other- 
wise; and  Wj  equals  | zPj  — A \{pj)z  — \2(pj)zS2(pp  — \3(pj)zS3(Pj)\  if  this  last  number  is  greater 
than  wcut , and  wcut  otherwise.  We  notice  that  given  j,  1 < j < m,  if  Wj  is  greater  than  wcut 
and  zs  is  set  to  z then  the  jth  term  of  Function  7.1  becomes  ( Wj/ujj)w 2 which  reduces  to  WjWj , 
i.  e.  Wj\zP]  — X\ (Pj)z  — \2(pj)zS2(Pj)  — X3(pj)zS3(Pj)\. 

The  derivative  of  Function  7.1  with  respect  to  zs  is 

m 

2EK'M')fe'  - Xl(pj)zs  - X2{Pj)zS2(pp  - A3 (pj ) zS3 (pj ) ) ( — Ai (pj ) ) . 

3= 1 

Setting  this  derivative  to  zero  we  are  then  able  to  solve  for  zs  and  thus  obtain  the  value  for 
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at  which  Function  7.1  is  minimized: 

m m 

2s  = (E(^M')(^  - A2 (Pj)z82(Pj)  - HPj)Zs3(p]))MPj))/Y,(Wj/{bj)MPj))2- 

j= 1 J=1 

The  following  iterative  procedure,  which  is  based  on  the  above  formula,  can  be  used 
for  attempting  to  minimize  approximately  Function  3.1.  Besides  STARGRIDPOINTS  and 
BARYCENTRIC,  the  procedure,  called  L1FIT3,  requires  a primitive  called  MAXRESIDUAL. 
It  computes  wcut  which  is  the  largest  absolute  value  of  a residual , i.  e.  the  largest  absolute 
value  of  a number  of  the  form  Z(p)  — \\\  ■ Z(Qi)  — W2  ■ Z(Q^)  — W3  • Z(Q3),  where  p is  in  /?,  Qi, 
Q-2 , Q 3 are  the  vertices  of  a triangle  in  T that  contains  p , and  W\ , W2,  W3  are  the  barycentric 
coordinates  of  p relative  to  the  triangle  so  that  p equals  W\  • Q\  + W2  • Q2  + FF3  • Q%.  The 
variable  wcut , which  we  call  the  residual  tolerance , is  used  throughout  the  execution  of  the 
procedure  for  the  purpose  of  avoiding  or  postponing  division  by  absolute  values  of  residuals 
that  are  less  than  its  current  value.  This  value  is  progressively  reduced  by  dividing  wcut 
before  each  pass  through  the  vertices  by  a variable  called  wscl.  This  last  variable  must  be  set 
on  input  to  a number  that  is  greater  than  1.0  and  preferably  less  than  or  equal  to  2.0.  The 
outline  of  the  procedure  follows. 

procedure  LlFIT3(i?,  A,  5",  T,  F,  Z,  Z,  n,  tol , npmax , u,  v,  wscl) 

begin 

flag  :=  0;  npass  :=  0; 
wcut:=MAXRESIDUAL(F,  Z,  Z); 

while  ( flag  = 0 and  npass  < npmax ) do 

begin 

flag  :=  1;  npass  :=  npass  + 1; 

wcut  :=  wcut/ wscl; 

if  ( wcut  < tol ) then  wcut  :=  tol ; 

for  z 1 until  n do 
begin 

if  (F(i)  G S \ S')  then 
begin 

(G,  m):=STARGRIDPOINTS(F(z),  R,  T); 
zold  :=  Z(F(i)); 

numerator  :=  0.0;  denominator  :=  0.0; 

for  j 1 until  m do 
begin 

(02,03,  Wuw2,  W3):=BARYCENTRIC(T.  G(j),  F(i )); 

IV  :=  Z(G(j))  - W,  • zold  - W2  ■ Z(Q2)  - W3  ■ Z(Q3 ); 
if  {W  > 0.0)  then  W :=  w else  IT  :=  n; 
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W :=  \W\; 

if  (W  < went)  then  W :=  went, ; 
numerator  :=  numerator + 

(wywo  • (Z(G(i))  - • Z(Q2)  - wy  Z(Q3))  ■ Wi; 

denominator  :=  denominator  + {W /W)  ■ (W\)2 

end 

2:new;  :=  numerator / denominator, 
if  (|  20/d  — znew  \>  tol ) then 
begin 

flag  :=  0; 

Z(F(i))  znew 

endendendendend 
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